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To locate the position and characterize the dynamics of a vacancy in a crystal, we 
propose to represent it by the ground state density of a quantum probe quasi-particle 
for the Hamiltonian associated to the potential energy field generated by the atoms 
in the sample. In this description, the /i 2 /2/x coefficient of the kinetic energy term 
is a tunable parameter controlling the density localization in the regions of relevant 
minima of the potential energy field. Based on this description, we derive a set 
of collective variables that we use in rare event simulations to identify some of the 
vacancy diffusion paths in a 2D crystal. Our simulations reveal, in addition to the 
simple and expected nearest neighbor hopping path, a collective migration mechanism 
of the vacancy. This mechanism involves several lattice sites and produces a long 
range migration of the vacancy. Finally, we also observed a vacancy induced crystal 
reorientation process. 
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I. INTRODUCTION 



Far from the melting, diffusion in solids arise when atoms migrate into empty sites in 
the crystal, leaving other empty sites behind them in which atoms can again migrate. Each 
of these empty sites is referred to as a vacancy, and the corresponding diffusion mecha- 
nism is said vacancy-driven. The atom/vacancy diffusion is a thermally activated process 
requiring the system to overcome free energy barriers separating the initial from the final 
state. In most cases, these free energy barriers are much higher than the thermal energy, 
therefore this process is a rare event, i.e. an event occurring with a frequency that is too 
low to be sampled by "brute force" Molecular Dynamics (MD) or Monte Carlo (MC) sim- 
ulations. Its study requires special simulation techniques, such as Temperature Accelerated 
MD(TAMD)/Temperature Accelerated MC (TAMC),™ MetadynamicP 31 or Adiabatic Free 
Energy Dynamics (AFED) 5 to efficiently explore the free energy surface of the system, or 
the string method in collective variables 6 to identify statistically relevant paths. These 
approaches require suitable collective variables (CVs) to describe the rare event. 

If one is not interested in the free energy, other methods exist that do not require CVs, or 
for which the identification of good CVs is not crucial. These methods can be classified in two 
groups: methods to identify the transition mechanism in the configuration space (e.g. the 
Nudged Elastic BandJ^the dimer method,^ the Transition Path Samplingp^), and methods 
for exploring potential energy surfaces, searching for their mechanical equilibrium points and 
saddle points (e.g. the Temperature Accelerated Dynamics, 11 the Hyperdynamics,^ and the 
ART met hod I ^IMJ) However, also these methods suffer from limitations (in addition to that 
already mentioned of not allowing to reconstruct the free energy). As for the first group, they 
require the a priori knowledge of the initial and final states of the process, and a tentative 
path belonging to the "reactive channel" one is interested in. As for the second group, 
despite significant progresses) 13 1 14 1 typically their reliability degrades with the complexity of 
the system. It must be stressed that these methods have been widely used to study defects 
(vacancy, interstitial, etc.) migration processes JSHID especially in "simple" systems. It must 
be also stressed that having a simplified but more expressive description of the vacancy 
migration process, like that attainable by CVs, is desirable as it allows to provide a more 
thoroughly physical picture of the events of this type, for example by effectively representing 
cooperative or collective effects. 
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The identification of good CVs for modeling the vacancy migration is still only roughly 
solved even for the most simple mechanism: the "local" vacancy migration mechanism. This 
mechanism consists in the hopping of the vacancy into a nearest neighbor lattice site. In 
this simple case, CVs have been proposed, e.g. in Refs. [TBI EH1 and [20j but turned out to 
be not completely satisfactory. 20 Moreover, more complex mechanisms have been suggested 
to take place in Bravais lattices and lattices with basis. For example, Da Fano and Jacucci 
found that the modeling of high temperature diffusion in Al and Na requires the inclusion of 
the vacancy double jump mechanism p2 Another case of complex vacancy migration mecha- 
nism is the collective migration of proton vacancies in hydrogen bonded network materials. 
For example, it was found that the collective proton transfer in [dabcoHj + fReOJ - (dabco 
= diazabicyclo[2.2.2]octane) is at the basis of the formation of the ferroelectric phase of 
this material.^ Vacancy migration in this material is coupled with the dynamics of the 
molecules forming the hydrogen bond network. Thus, a suitable CV must be able to take 
into account this phenomenon. In even more complex crystals, such as clathrate hydrates of 
molecular gases, the vacancy is associated to (missing) guest molecules (CH 4 , CO2, etc.) and 
their diffusion mechanism requires the cooperative effect of the water molecules forming the 
framework of the crystal.^ Also in this case, the CV should be able to take into account the 
concerted dynamics of the framework and guest molecules. In materials with a high concen- 
tration of vacancies, such as yttria and scandia stabilized cubic zirconia, vacancy- vacancy 
correlations were found to play a crucial role in the mass transport mechanism (for example 
of oxygen in the case of zirconiaP^and using methods based on the configuration space it 
might result difficult to characterize this phenomenon.^ Finally, the vacancy diffusion is 
proposed to be the limiting step in the nucleation and growth of some nanostructures, for 
example in CU2S nanowiresp3 Studying the mechanism of this process requires a description 
of the vacancy, and its migration, that is independent on the orientation of the growing nu- 
cleus (which is not known in advance), robust with respect to some degree of disorder that 
might be present in small crystal-like nuclei, or even to the change of crystalline symmetry 
along its growth. 

While the description of the simple local migration mechanism in Bravais lattices can be 
obtained by some ad hoc CVs, the treatment of the more complex cases mentioned above 
requires the introduction of "general" CVs to represent vacancies and their dynamics. By 
"general" we mean CVs that are not specific to any crystal symmetry or orientation, or 
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tailored to monitor a specific migration path. Another requirement is that the CVs work 
also with "stressed" crystalline systems, such as nano crystals, that might by non-uniform 
or non isotropic (see, for example, Refs. I2TH291) . The aim of this paper is to introduce such 
a general set of CVs. We illustrate the use of our CVs by investigating vacancy migration 
paths in a 2D crystal of purely repulsive Lennard- Jones particles (Week-Chandler- Andersen 
particles^ - WC A) . We believe that our CVs can be used also in more complex systems of 
the kind mentioned above. 

The remainder of this paper is organized as follows. In Sec. [II] we introduce the general 
formalism of our description of the vacancy, and derive a set of CVs that are numerically 



more convenient to study the vacancy migration process. In Sec.yTywe test our approach by 
first characterizing the vacancy migration observed during a high temperature MD trajectory 
of a simple 2D WCA crystal, and then using a reduced CV set in combination with TAMC 
to find interesting migration paths (and other events) in the same system. In Sec. IV we 
draw some conclusions. 



II. THE QUANTUM PROBE QUASI-PARTICLE REPRESENTATION OF 
A VACANCY. 

In this section we introduce the observable that we consider the most adequate to rep- 
resent a vacancy and track its dynamics. We postpone to Appendix [A] the description of 
simpler alternatives that, although more intuitive to grasp, are not completely satisfactory. 

A vacancy is the lack of an atom at a lattice site in a crystalline system. Therefore, 
describing a vacancy amounts to identifying the empty site in a crystal and following its 
evolution in time. Our starting point is to represent a vacancy as a quantum probe quasi- 
particle subject to the field produced by the atoms in the sample: 

V(x;R) = f>(|R,-x|) (1) 

1=1 

where x is the position of the probe particle, R is the 3N vector of the atomic positions 
(Hi is the position of the atom i) and v(r) is the pair potential governing the interaction of 
the atoms in the system. The extension to more complex interaction potentials, including 
ab initio ones, is conceivably feasible (see Appendix [B]). The potential energy in Eq. [j] 
is that of a probe particle located at x when the atoms are in the configuration R. To 
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FIG. 1. Series of atomic configurations, and the corresponding V(x;R) potential, along the local 
migration path in a 2D trigonal lattice of WCA particles (see Sec. |III A| for all the details). F(x; R) 
is reported in Lennard- Jones units (these units are used throughout the text) and represented 
according to the colormap reported in the bar at the bottom of the figure. The simulation box 
is also reported in the figure (periodic boundary conditions are used in all calculations discussed 
in the text). The top panel corresponds to the configuration in which the vacancy is located on 
a lattice site. The arrow indicates the path followed by the moving atom during the (ideal) local 
vacancy migration process described in the text. In the next two panels are shown the atomic 
configuration and the potential V (x; R) at intermediate states. In these states, the potential 
V(x; R) presents already the characteristic "double- well" shape. In the bottom panel is reported 
the mid- way configuration and the corresponding potential V^(x; R). In this state, the two wells of 
V(x; R) are symmetric. 
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give to the reader an intuitive argument of why the above potential is a key element in 
the vacancy description, in Fig. [j] we report the V(x;R) along the (ideal) local migration 
process in a 2D WCA trigonal crystal (more details on the calculation of V^x; R) are given 



in Sec. III). In this process one atom, nearest neighbor of the vacancy, moves along the linear 
path connecting its (initial) lattice and the vacancy site. At the beginning (top panel), the 
potential is characterized by one deep minimum in correspondence of the vacancy, and other 
local minima located in between the atoms. From now on, we will refer to the former as 
vacancy minimum (minima) and to the latter as crystal minima, as they are present also 
in perfect crystals.^ While the atom moves toward the vacancy (the two central panels), a 
second vacancy minimum is formed at its original site. When the atom is mid- way along 
the path (bottom panel), the potential V(x; R) presents two vacancy minima of equal value 
in the region of the two lattice sites involved in the process. 

One could think of using directly the potential V(x; R) to describe the vacancy. However, 
this representation is way too complex, since V^x; R) is typically characterized by very many 
minima of which only one or few of them are relevant for the description of the vacancy 
migration process. We propose to describe the vacancy by the ground state probability 
density of the quantum probe quasi-particle, p(x|R) = |?/;(x;R)| 2 , where ^(x;R) is the 
function minimizing the Ritz functional 

<Kx- R) - ar £ min W*;R)|K(x;RMx;R)) (2) 

In Eq. (2] 7{(x;R) = — aV 2 + V(x;R), and a > is a tunable parameter we will dis- 
cuss shortly.^ ^(x; R) can, equivalently, be defined as the ground state of the Schrodinger 
equation associated to the Hamiltomian H(x;R), but the formulation in terms of the Ritz 
functional makes more simple the analysis of some of its properties. 

In the following we will show that, for a suitable value of a (discussed in detail in 
Sec. HI Ap , p(x|R) is localized around the vacancy minima of the V^xjR) potential and 



is smooth. To understand why p(x|R) is localized around vacancy minima, consider sepa- 
rately the potential and kinetic energy contributions to the expectation value of the operator 
H(x; R). The potential energy part can be recast into the form (^(x; R)|V(x; R)|^(x; R)) = 
/ dxl/(x; R)p(x|R). This term is minimized by a density, and the corresponding wavefunc- 
tion, that is localized around the minimum of V (x; R) (p(x|R) = 5(x — x(R))). However, 
to a very localized wavefunction corresponds a very high "kinetic energy" , as can be seen by 



noticing that the "uncertainty principle" imposes that — (^(x; R)|aV 2 ^(x; R)) > a/4(Ax 2 ), 
where (Ax 2 ) = (x 2 ) — (x) 2 P^ Therefore, the wavefunction minimizing the Ritz functional 
above is the trade off between the need to be localized around the minimum (minima) of 
the potential V (x; R) and the need to be not too localized as, otherwise, the kinetic energy 
would be too large. The parameter controlling the degree of localization of the wavefunction 
is a. If a is small, V^x; R) is the dominant term and the ground state density is very local- 
ized around the absolute minimum of the potential. If, on the contrary, a is large the kinetic 
energy becomes the dominant term and the ground state density must be highly spread out. 

Typically, the kinetic energy of the ground state has to be not too large. This fact has 
a crucial consequence on the ability of the quantum density representation to distinguish 
between vacancy and crystal minima also in those cases in which the values of the potential 
in the two regions are not well separated. To show this, we start by noticing that crystal 
minima are narrower than vacancy minima. Therefore, the wavefunction corresponding 
to a p(x|R) that presents maxima at crystal minima must grow from zero (at the atomic 
positions) to a large value and go back to zero on a small length scale, corresponding to 
the inter-atomic distance. A wavefunction of this type will have a large gradient in that 
region and, therefore, a large kinetic energy, and will not correspond to the ground state of 
the above Hamiltonian. Therefore, the ground state density will be peaked at the (wider) 
vacancy minima, also when the value of V (x; R) in the two different types of minima is not 



very different. Indeed, in Sec. Ill we will show that p(x|R) is localized at the vacancy even 
in the case in which the values of the potential at crystalline and vacancy minima are the 
same. 

The kinetic energy term is also responsible for the smoothness of the ground state density. 
In fact, using the same argument used for its localization, to a rough wavefunction would 
correspond a high kinetic energy. Thus, the ground state wavefunction must be smooth and, 
therefore, the ground state density will be smooth too. 

Summarizing, for a suitable choice of a, the quantum density is localized around the 
vacancy minima of the potential ]/(x;R), even when the value of the potential at the 
vacancy and crystal minima is similar, and smooth. 

The discussion above, and the computer experiments reported in Sec. |III A[ indicate 



that p(x|R) is able to localize the vacancy and to characterize its migration path. There- 
fore, we could think of using p(x|R), or better its value on a discretization of the x-space 



({p(xi|R)}i = i ? M)j as a vectorial (field-like) CV in rare event simulations. However, this 
approach presents problems. First of all, the dimensionality of this vectorial CV, which 
corresponds to the number of grid points in the discretization of the x-space, would be quite 
large, typically much larger than the number of atoms in the simulation sample. This is 
because we need to characterize accurately the density between atoms. For example, in the 
simple local migration event shown in Fig. [I] we need to represent accurately the p(x|R) in 
the region between the moving atom and its nearest neighbors, where a new minimum of 
the potential V (x; R) is forming. This requires to have a mesh with several grid points per 
lattice site, from which we can conclude that the typical dimensionality of the discretized 
field-like CV proposed above is larger than the number of atoms. It is also worth mentioning 
that, thanks to the localized nature of the p(x|R), only a small subset of the elements of 
this vectorial CV is non negligible at each atomic configuration. In other words, the amount 
of information provided by {p(x^; R)}i=i,M can be redundant if we just want to identify the 
"position" of the vacancy and follow its dynamics. 

The use of {p(x^; R)}^=i 5 m poses also another, severe, problem: the functional set of 
the ground state densities of the Hamiltonian / H(x;R) is, in general, unknown. This prob- 
lem, referred to as a ]/-representability of the ground state densities" in density functional 
theory,^ adds an additional difficulty in using this CV in both guided (umbrella sampling, 35 
blue moon, m etc.) and unguided (TAMD/TAMCF 2 Metadynamics, 3 4 etc.) rare event 
simulations. In the first case, because we do not know how to set l/-representable re- 
straint/constraint values of the CV. In the second case, since the random value of the CV 
generated along the biased dynamics might be non l/-representable, the simulation will spend 
a non negligible amount of time sampling values of {p(x^; R)}^=i 5 m that are not in the set 
of the possible solutions, thus reducing the efficiency of the approach. A possible, still infor- 
mative, simplification of the representation consists in using as CVs few low order moments 
k u of the density p(x|R), where n is the order of the moment. While the full representa- 
tion of the density requires the infinite set of its moments {^ n } n =i,oo^ 7 its low order terms, 
/^(R) = x(R) = fdx. p(x|R)x and ^(R) = c^(R) = J dx. p(x|R)(x a - x a )(x p - x p ), 
with a,/? = lj...,d(d dimensionality of the space), seem adequate to catch the features of 
p(x|R) characterizing the state of the vacancy. The first moment x(R) tells us where the 
vacancy is located, whether on a lattice site or in between sites, the latter case indicating 
that a migration event is taking place. Taken alone, x(R) is not sufficient to characterize 
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the vacancy migration process. For example, the x(R) cannot distinguish between a local 
vacancy migration event, involving only one atom, from a collective one, with several atoms 
moving at the same time from their initial site to the next one, resulting in a long range 
migration of the vacancy. The second moment can help distinguishing between these two 
kinds of processes because c 2 (R) measures the width of the vacancy, which depends on the 
number and positions of the sites involved in the process. In particular, the trace of this 
matrix, Tr[c 2 (R)], is small if the migration is local and large if the migration is long range. 

In conclusion, the collective variables we propose to use to study the vacancy migration 
process are the first moment and the trace of the second moment matrix of the density 
p(x|R). Should this set of collective variables result insufficient, e.g. should the anisotropy 
or the asymmetry of the p(x|R) be a feature characterizing the vacancy migration in the 
system at hand, the set could be improved by replacing the trace of the second moment with 
the entire matrix (or its eigenvalues), or by adding moments of higher order, for example 
the skewness. The full approach, i.e. the use of p(x|R) as CV, would still be desirable, but 
this requires further, nontrivial, thinking to solve the problems mentioned before. 



III. RESULTS AND DISCUSSION 



This section is divided in two subsections. In Sec. |III A| we show that the quantum 
density representation is able to describe the vacancy migration process at zero and finite 
temperature, both when the potential in the crystal and vacancy minima has well separated 
values (the WCA potential will be used as a representative example) and when these values 



are superimposed (hard disks). In Sec. Ill B we study the vacancy diffusion in a 2D crystal 



using the reduced set of collective variables introduced at the end of Sec. [TTJ and exploit 
these data to identify possible diffusion paths. 

A. Use of the quantum probe quasi-particle representation for the 
identification of a vacancy in a 2D crystal. 



The WCA pair potential is a purely repulsive potential of the form: 



v(r) 



4e 



(a/r) 12 - (a/rY 



€ Vr < 2 1 / 6 <7 



0, Vr > 2 1 /V 



(3) 



For uniformity of notation, we will denote by a also the size of the hard disks. The simu- 
lations are performed on a sample of 24 atoms and one vacancy in a periodic trigonal 2D 
box, corresponding to one defected 5x5 layer of the (111) surface of a face centered cubic 
lattice (see Fig. [j]). For the WCA system, the lattice constant was fixed to 1.075a, to be 
compared with 2 1 / 6 a ~ 1.222 a, the value at which the WCA potential becomes zero. For 
the hard disk system, we set the lattice constant to 2.5 a. 

Before moving to the description of the calculation of p(x|R) we must explain how to set 
the value of a. In Sec. |TT] we explained that a must be neither too small, as otherwise the 
p(x|R) will be too rough and could not distinguish between vacancy and crystal minima, 
nor too large, as otherwise the density will tend to be uniformly distributed over the entire 
x-space. In practice, we set its value such that the width of p(x|R) (Ax = yj ((x — (x)) 2 )) 
in a crystal with a vacancy at the equilibrium configuration is of the order of the typical 
interactomic distance in the crystal, e. g. the size of the unit cell for Bravais lattices, or the 
Van der Waals radius of the missing atom for more complex crystals. In the present work 
we set a = 1500. Below, we investigate the effect of a on the p(x|R) and show that this 
value is adequate for the WCA system. 

The ^(x; R) is computed by expanding the wavefunction on a plane wave basis set of ele- 
ments Xi(x) = 1/y/V exp[iG^x] (see Ref. [38] for technical details of planewave calculations). 
The wavevectors G^ satisfy the usual condition that the corresponding Xi( x ) is periodic over 
the simulations box, which amounts to set the {Gi}^ to the points of the reciprocal lattice of 
the simulation box (G$ = 27iH~ 1 \i i) where H — [a b], with a and b representing the edges 
of the simulation box, and ui = ±1, ±2, . . .). The expansion is limited to the planewaves 
satisfying the condition a|G^| 2 /2 < E cut . E cut is fixed such that the ground state eigenvalue 
of the Hamiltonian matrix at several atomic configurations is well converged, i.e. the varia- 
tion with E cut is less than 10 -3 Ae, with Ae = ei — eo and eo and e\ the eigenvalues of the 
ground and first exited state of the Hamiltonian, respectively. In the present case, E cu f was 
set to 5 x 10 6 , corresponding to a 100 x 78 points Fourier mesh. 

The ground state of H(x; R) in the T-point approximation is obtained using the Lanczos 
iterative method.^ In particular, we used the implicitly restarted version of the met ho d,^^ 
which is more efficient than the original one. 

The fact that the WCA and hard disk potentials diverge at x = R^ and |x — R^| < a, 
respectively, poses a problem for the calculation of the matrix elements of the Hamiltonian 
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in the planewave basis set. This problem was solved by replacing the original pair potentials 
with non diverging approximations. In the WCA case, the original pair potential is replaced 
with the following approximation: 



We tested the dependency of our results on r cut and verified that for a small enough value 



independent from the chosen value of r cut (in the present work r cut was set to 0.35). If a 
more accurate approximation were needed, it would be possible to develop a pseudopotential- 
like approach.^ For the hard disk case, we solved the problem by setting the value of the 
potential to a large but finite value for r < a. In the calculations discussed in the following, 
we set this value to 10. Also in this case, we tested that above a given threshold our results 
are independent on the specific value chosen. 

In our calculations we also need to ensure that the vacancy wavefunction does not inter- 
acts with its periodic images. To test this, we computed the dependency of p(x|R) on the 
sample size for selected configurations along the ideal local migration path (Fig.[T]). We mea- 
sured this dependency by computing two quantities: i) max x Ap z>a (x|R)/ max x p 5x5 (x|R), 
where Ap lxl (x|R) = |p z>a (x|R) — p 5x5 (x|R)| and p 2XZ (x|R) is the vacancy density in a i x i 
system, and ii) / rfx Ap zxz (x|R). For i = 10 — 30 (with increments of 5) we found that 
max x Ap ixi (x|R)/max x p 5x5 (x|R) < 5xl0~ 4 and / rfx Ap ix '(x|R) < 10" 5 . The very small 
dependence of p(x|R) on the size of the simulation box is due to its very localized nature. 

We start the presentation of our results by showing that the quantum density description 
is able to localize the vacancy in the correct region of the space for the well ordered atomic 
configurations reported in Fig. [j] (corresponding to T = 0), both in the case of the WCA 
and hard disk pair potentials. In Fig. [2] we plot the p(x|R) associated to the configurations 
of Fig. [IJ Initially, when the atoms are at their equilibrium position (top panel), the density 
is localized at the vacancy site. Then, while the migrating atom proceeds along its linear 
path, the density follows the vacancy minimum of the V^xjR) (the two central panels), 
finally splitting in a symmetric bimodal p(x|R) when the migrating atom is in the mid- 
way configuration (bottom panel). This figure clearly shows that the quantum density 
representation is able to describe the vacancy all along the (ideal) local migration path in a 





of this parameter, such that V^x; R) = V^(x; R) around the vacancy minima, the results are 



system in which there is a sizeable difference in the value of V (x; R) between vacancy and 
crystal minima. However, as explained in Sec.[TTJ the quantum density representation is able 
to characterize the vacancy also when there is no such a separation. This is shown in Fig. [3| 
where we report data analogous to those of Fig. [2] for a system of hard disks (only data for 
the equilibrium and mid- way configurations are shown). Also for the hard disk system the 
quantum density representation is able to properly describe the vacancy and its dynamics 
along the ideal local migration path. 

To be able to use the quantum probe quasi-particle description we need to show that it 
is possible to find a values of the parameter a such that p(x|R), and its low order moments 
considered before, are able to identify the vacancy position and follow its migration. This 
is verified using x(R) and Tr[c 2 (R)] computed at selected configurations along the local 
migration path (Fig. [I]). The criterium we adopt is that the description is valid if the 
variation of the values of these observables computed at significantly different configurations 
along the migration paths (e.g. the initial and the mid- way configuration of Fig.[T]) are larger 
than thermal fluctuations. For the WCA system at T = 0.5, for example, the root mean 
square thermal fluctuation of the first and second moment, Ax(R) and ATr[c 2 (R)], are 0.1 
and 0.05, respectively (see left panels of Fig. [8]). In Fig. [I] we report x(R) and Tr[c 2 (R)] vs 
a. A curve in the top, central and bottom panel of this figure represents xi(R), x 2 (R) and 
Tr[c 2 (R)], respectively, at a given atomic configuration vs a. Curves with same color and 
symbol in different panels refer to the same configuration. We identify three zones, denoted 
1, 2 and 3 in the figure. In zone 1 p(x|R) is essentially independent on the value of a, and 
x(R) and Tr[c 2 (R)] are almost constant over the entire zone. This happens when the p(x|R) 
is more localized than the interatomic distance (1.075 a in our case) in the vacancy minima. 
In this zone, the difference of the values of x(R) and Tr[c 2 (R)] between the initial and 
mid- way configuration are larger than their thermal fluctuation at T = 0.5. Moreover, one 
of the two components of the first moment, ^(R), is well separated also for configurations 
closer to each other, e.g. at 1/5 and 2/5 of the arclength distance I along the local migration 
path, which is enough to distinguish between these two different states. In zone 2, p(x|R) 
is no longer independent on a and, for a given configuration, x(R) and Tr[c 2 (R)] change 
continuously with a. However, also in this case, for a given value of a, the values of x(R) 
and Tr[c 2 (R)] at significantly different configurations along the local migration path are well 
separated. Finally, in the zone 3, p(x|R) is no longer localized in the vacancy region and, 
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FIG. 2. V(x;R) (left) and p(x|R) (right) for four atomic configurations along the ideal local 
vacancy migration path in a system of WCA particles. 

as a consequence, x(R) and Tr[c 2 (R)] take values essentially independent from the atomic 
configuration. Obviously, values of a belonging to this range are unsuitable. Summarizing, 
the quantum probe particle description is suitable to represent the vacancy and track its 
dynamics over a wide range of values of a (i.e. the choice of the value of a is not critical), 
and the criterion we mentioned before, that a must be such that Tr[c 2 (R)] x / 2 is close to the 
interatomic distance, is adequate. 

In order to test whether this model is still able to represent the vacancy and its dynamics 
when the system is at finite temperature, we run a high temperature (T = 2.5) MD simu- 
lation with a Langevin thermostat. In this case, we observed one vacancy migration event 
in a 10000-step long simulation. On the short timescale of our MD simulation the system 
remains crystalline, even though its equilibrium state at this T could be the liquid one. In 
the left column of Fig. [5] is reported the potential V^x; R) along the (local) migration event 
mentioned above. By comparing the potential shown in Fig. [5] and Fig. [I] we see that, at 
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FIG. 3. F(x; R) (left) and p(x|R) (right) for a system of hard disks in the equilibrium and mid-way 
configurations along the ideal local migration path. The blue color in the left-hand panel denotes 
V (x; R) = 0; the red V (x; R) = oc. 

finite temperature, the V (x; R) along a local migration event has still a shape similar to 
that at T = 0. However, few narrow minima, of magnitude similar to those of the vacancy 
ones, are also present. These minima are due to large displacements of atoms not involved in 
the vacancy migration process out of their equilibrium position. Despite the more complex 
shape of the potential in the configurations visited at finite (and high) temperature, the 
p(x|R) is still able to correctly identify the position of the vacancy, both when the vacancy 
is located at a lattice site and when it is moving from one site to another (see the right-hand 
panel of Fig. [5]). Indeed, while some density leaks to the crystalline minima, the value there 
is negligible compared to the one around the vacancy ones. 

We also analyzed the high temperature MD trajectory to test the ability of x(R) and 
Tr[c 2 (R)] to monitor the vacancy migration processes, if our CVs are good, their oscillations 
when our system is in a metastable state should be much smaller than their variation when a 
vacancy migration takes place. This condition is indeed verified, as can be seen by comparing 
the typical oscillation of the CVs against their (simultaneous) change at the timestep ^ 4700 
in the x(R) and Tr[c 2 (R)] timelines along the T = 2.5 simulation (Fig. [6]), corresponding 
to the migration event. By visualizing the trajectory, we verified that this level is of local 
kind (this is, indeed, the event shown in Fig. [5]). In Fig. [6] we also notice several peaks in the 
Tr[c 2 (R)] curve to which do not correspond any significant change in the x(R). We analyzed 
the origin of this behavior by looking at the p(x|R) during one of these events, namely the 
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Zone 1 Zone 2 



Zone 3 



FIG. 4. x(R) and Tr[c 2 (R)] as a function of a at the configurations shown in Fig. [T] a is reported 
on a logarithmic scale. Vertical lines are inserted to help the comparison of values of x(R) and 
Tr[c 2 (R)] at a given value of a. I denotes the arclength distance along the ideal local migration 
path of Fig. [I] starting from the initial condition. Zone 1, 2 and 3 denote the three regions of a in 
which x(R) and Tr[c 2 (R)] have a different behavior (see text). 



one at the timestep ~ 5600 (indicated by an arrow in Fig. |6j). A series of snapshots taken 
along this event are shown in Fig. [7| The sudden change in the second moment is due to 
an unsuccessful migration event. This produces a broadening of the p(x|R), which partly 
populates also the minimum of the V (x; R) that is forming on the crystal site initially 
occupied by the migrating atom. However, the process does not end with a migration and 
the final p(x|R) is still localized on the initial empty crystal site. The fact that the first 
moment does not change (significantly) during this event is due to the combination of two 
factors. On the one hand, during the attempted migration event there is an increase of 
the density on the site of the moving atom. This should move the center of the density 
p(x|R), i.e. x(R), toward this site. However, at the same time, the maximum of the density 
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FIG. 5. F(x; R) (left) and p(x|R) (right) at few selected atomic configurations along the migration 
event observed during an unbiased simulation at T = 2.5. The arrows in the left column point to 
representative narrow minima of the potential y(x;R) not associated to the vacancy (see text). 



on the original vacancy site moves in the opposite direction, following the position of the 
original vacancy minimum of the potential V(x;R), which is "pushed" in this direction by 
the moving atom. This second effect would move the x(R) in the opposite direction and the 
two effects, essentially, compensate. This can be seen in Fig. [7| where, together with the 
snapshots of the potential V^x; R) (left) and the density p(x|R) (right) along the attempted 
migration event, we report the first moment x(R) (denoted by the white cross), see also the 
discussion in the caption. 

In conclusion, the high temperature MD test indicates that x(R) and Tr[c 2 (R)] are both 
needed to monitor the vacancy migration process. At the same time, this test makes us 
confident that these three CVs are sufficient to study this process using rare event techniques. 
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FIG. 6. The two components, x\ and #2, of the first moment (top) and the trace of the sec- 
ond moment (bottom) of the density p(x|R) as measured along a high temperature (T = 2.5) 
simulation. 



B. Vacancy diffusion path in a 2D WCA crystal by Temperature Accelerated 
Monte Carlo. 



In this section we report the results of rare event simulations aimed at identifying vacancy 
migration paths in a 2D WCA crystal at finite temperature using x(R) and Tr[c 2 (R)] 
as CVs. These calculation will be performed using TAMC.™ In TAMC the dynamical 
system consists of the original atomistic variables plus a set of extra variables z = {^}i=i, m , 
associated to a set of collective variables {0i(R)}i= m , x(R) and Tr[c 2 (R)] in the present case. 
The two sets are coupled via the potential energy term C4(z, R) = J2i=i ™ ^/2(£^(R) — Zi) 2 . 
The atoms and the z are "evolved" together, the atoms according to a standard Metropolis 
Monte Carlo governed by the physical potential energy plus the term C4(z,R), and the 
z according to a constant temperature dynamics (Langevin dynamics in the present case) 
governed by the potential C4(z, R). The inertia of the z, a tunable parameter in this method, 
can be set such that the evolution of these variables is adiabatically separated from the Monte 
Carlo on the atoms, so that the latter samples the conditional probability density function 
m(R|z). It can be shown (see Ref. 2) that under these conditions the variables z sample the 
probability density function P$(z[) = exp[— /3*F&(z)], where /3* = (Jcb Boltzmann 

constant) and, for k large enough, F&(z) is the free energy at the physical temperature T 
associated to the state 0(R) = z. T* is the temperature of the z variables, which can be 
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FIG. 7. Potential F(x;R) (left) and density p(x|R) (right) computed on configurations corre- 
sponding to few snapshots along the avoided migration event at the timestep ~ 5600 of the high 
temperature MD simulation. The white cross indicates the position of x(R). When the vacancy is 
located at a crystal site (top and bottom panels), x(R) is approximately in correspondence of the 
maximum of the p(x|R) and the vacancy minimum of the potential T^(x; R). During the attempted 
migration event (central panels), the main peak of the density is shifted toward the bottom-left 
with respect to x(R), while a long tail is formed at its top-right. 

different from the physical temperature. By defining T* such that the associated thermal 
energy is higher than the barriers separating free energy minima, we are able to efficiently 
explore the free energy surface. 

In this work, we want to identify possible vacancy migration paths in a 2D WCA crystal 
by direct inspection of microscopic configurations explored in our TAMC runs. The recon- 
struction of the free energy surface, the accurate determination of the migration paths, and 

18 



the corresponding rates via the transition state theory with dynamical corrections) 11 will 
be the objective of a forthcoming study on a more realistic 3D system. 

We performed a 2.2 x 10 6 steps TAMC simulation on the already described 2D WCA 
crystal with one vacancy. The physical temperature was set to T = 0.5 and the temperature 
of the CVs to T* = 8.5. The simulation was started from a configuration in which the 
vacancy and all the atoms are at the lattice sites. In Fig. [8] we compare the values of x(R) 
and Tr[c 2 (R)] as obtained from a MC simulation at T = 0.5 with those obtained from the 
TAMC simulation. We did not observe any migration event in the MC simulation. This 
is consistent with recent result^ on vacancy migration in a closely related 2D crystal at 
T = 0.1, for which the energy barrier of the local migration process was estimated to be ^ 20. 
At variance with MC simulations, in the TAMC run we observed several migration events, 
identified by the simultaneous change of x(R) and Tr[c 2 (R)]. In the TAMC simulation, 
we also observe many peaks in the Tr[c 2 (R)] not associated to any change of x(R). As 
explained in the previous section, these peaks are due to aborted migration events, which at 
T* = 8.5 are quite frequent. By visual inspection of the atomic "trajectory" of the TAMC 
simulation, we identified three different kinds of events described in the following. The first 
one corresponds to a local migration event. The atomic trajectory along one of the eight 
events of this type observed in our simulations, namely the one at the step ^ 1.1 x 10 6 (see 
Fig. [8]), is shown in Fig. [9]/ A. It can be noticed that while the atom 1 is migrating toward 
the vacancy the next neighbor atom (atom 2 in the figure) is moving in the same direction, 
giving rise to a concerted motion. This kind of synchronous motion is found also in the 
other local migration events observed in the TAMC simulation. 

In Fig. [9J/B we report another kind of event, corresponding to the step - 2 x 10 6 of the 
simulation, that we call non local migration event. In this case, several atoms and lattice 
sites are involved in the process, namely four atoms and five lattice sites. The trajectory 
reported in the figure shows that this event does not consist of a series of independent local 
migration events close in time. Rather, the atoms move all together and the process ends 
when they have all reached the final lattice position. The result of this event is that the 
vacancy is transferred at a distance of four lattice sites from its initial position. This is, 
indeed, a concerted "multiple jump" process of the kind identified by Da Fano and Jacucci 
in high temperature Na and Al samples.^ During the TAMC simulation, we observed only 
one event of this type, against eight local migration events. This is consistent with the 
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FIG. 8. Comparison of the first moment and the trace of the second moment of the density p(x|R) 
as obtained from (unbiased) MC at T = 0.5 (left) and TAMC at the same T and at T* = 8.5 
(right). The first component of the x {x\) is denoted by a continuous - black - line, while the 
second component (^2) by a dotted - red - line. The dashed ellipses indicate events that are 
discussed in the text. The inset shows a zoom of the Tr[c 2 (R)] along the TAMC simulation in the 
interval in which the system was in a reoriented crystalline state (see text). The red line correspond 
to the bottom of the Tr[c 2 (R)] curve in the rest of the simulation. 
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FIG. 9. Trajectories along a local (A) and non-local (B) migration event. The particles are colored 
in red at the beginning of the trajectory and, passing by green, become blue at the end. The 
boundary of the simulation box is also reported. 
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FIG. 10. Atomic configurations corresponding to the misoriented crystal. The stress due to wrong 
orientation is accommodated by the formation of a wide vacancy defect. During the interval 
~ 2.3 — 5 x 10 5 TAMC steps the misoriented crystal kept the same geometry (trigonal) but its 
orientation and the position of the vacancy changed several times. The misorientation of the 
crystal is put in evidence by reporting on the plot the original lattice vectors a and 6, and the 
lattice vectors of the reoriented crystal a f and b' . 

observation of Da Fano and Jacucci that the non local migration process is statistically 
relevant only at high temperature (i.e. close to the melting). 

Finally, the third kind of process observed in our simulations corresponds to a reorien- 
tation of the lattice. This process, taking place at the step ~ 2.3 x 10 5 of the TAMC run, 
starts as a non-local migration event from a properly oriented crystal structure. This initial 
step is followed by a global change of the atomic positions into a lattice of analogous sym- 



metry (trigonal) but with a different orientation (see Fig. 10). This misoriented lattice is 



incommensurable with the simulation box, thus producing a large stress on the system. The 



system reacts to this stress by forming wider vacancies, such as those shown in Fig. [10| The 
formation of these defects is reflected on the value of the Tr[c 2 (R)]. In fact, in the inset of 
Fig. [8] you can see that in the interval ~ 2.3 — 5 x 10 5 the bottom of the Tr[c 2 (R)] curve is 
higher than in the rest of the simulation, when the system is in the ordinary orientation. The 
process ends with another global change of the atomic position restoring the original orienta- 
tion of the crystal. This process is most likely an artefact of the large disorder present in the 
crystal due to the high vacancy concentration in the sample, which is orders of magnitude 
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higher than the typical value in bulk systems. Nevertheless, it is very promising that our 
CVs are able to accelerate such a process. This could allow to investigate the formation of 
similar states/processes in low dimensional systems (e.g. interfaces, nanocrystals) in which 
the formation of misoriented crystals might be more favorable than in bulk systems. 



IV. CONCLUSIONS 



In this paper we have introduced a novel field-like observable able to locate a single va- 
cancy in a crystal and follow its dynamics. This observable is the ground state probability 
density of a quantum probe quasi-particle for the Hamiltonian associated to the potential 
energy field generated by the atoms in the sample. To exploit this observable in practice, 
we derived from it a small set of collective variables that, used in conjunction with rare 
event techniques, allowed us to study possible vacancy migration paths in a 2D crystal 
of Week-Chandler-Andersen particles. Our simulations revealed, in addition to the sim- 
ple nearest neighbor vacancy migration mechanism, a long range migration path consisting 
of the simultaneous jump of several atoms. Moreover, we observed a crystal reorientation 
process induced by a multiple jump vacancy migration event. Work is in progress to gener- 
alize this description to the many vacancy case and to the possible creation/annihilation of 
vacancy /interstitial pairs. 
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Appendix A: Inadequacy of simpler representations. 



One might think of using simpler representation of the vacancy. In this appendix we 
describe two of them and discuss their inadequacy 

The first one consists in representing the vacancy as a classical probe particle. In this 
case, the vacancy is represented by a particle positioned at x(R) = argmin x V(x; R), i.e. 
the mechanical most stable equilibrium position of a probe particle in the field generated 
by the atoms. This works when the vacancy is at a lattice site but it is no longer adequate 
when a migration process takes place, as we illustrate by the following ideal experiment 
of a local vacancy migration event. For a 2D WCA crystals or similar systems, when a 
vacancy is located at a lattice site the corresponding minimum of V^x; R) is much deeper 
than the crystal minima (see top panel of Fig. [T]) and the classical probe particle description 
is univocal and adequate. When the atom is mid-way along this path (bottom panel of 
Fig. [lj/D), the potential l/(x;R) presents two minima of equal value in the region of the 
two lattice sites involved in the process. In this configuration, the definition of x(R) be- 
comes ambiguous. However, even before reaching this double-minimum state, the vacancy 
migration representation given by the classical probe particle description is unsatisfactory. 
In fact, while the atom moves along the path described above, a second minimum starts 
forming close to the original lattice of the migrating atom (two central panels of Fig. [I]). 
The magnitude of this minimum increases while the atom moves toward the empty site, 
until it matches the magnitude of the other minimum in the mid- way configuration, but 
x(R) remains unchanged. Thus, the probe particle description is unable to represent the 
continuous evolution from the initial to the mid- way configuration of the system. 

A natural improvement over the previous description could be obtained by describing the 
vacancy in terms of the (classical) probability density pc(x|R) to find a probe particle at x 
conditional to the atoms to be at the configuration R: 

n /„|td\ _ exp[-^(x;R)] 

PC(X|R) " /dsexph9V(x;R)] (A1) 



In Eq. Al f3 is a tunable inverse "temperature" controlling the localization of pc(x|R): the 
higher is /3 the more p(x|R) is localized around the deeper minimum/minima of the potential 
V(x; R). pc(x|R) is able to correctly represents configurations such as the one shown in the 
bottom panel of Fig. [TJ where the vacancy is split in two, when the difference between the 
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(A) (B) 

FIG. 11. 3D representation of the classical (A) and quantum (B) densities in finite temperature 
system at the same configuration of the third panel of Fig. [5j The classical density has been 
computed at f3 = 0.1, the lowest (3 with as associated density with negligible density on crystal 
minima. The scales of the panels (A) and (B) is set such that their maxima have a comparable 
magnitude. 



values of V (x; R) at vacancy and crystal minima is sizeable and when the system is at T = 0. 
However, when these two conditions are not met also the classical density representation 



has problems. To illustrate this, let us consider first the case discussed in Sec. |III A| of a 
system composed of hard disks. In this case the classical density representation would be 
completely inadequate since, due to the fact that the potential V (x; R) is zero everywhere 
apart within the disks, the density would be uniform outside the disks. 

Let us now consider finite temparature effects. In this case, atoms move out of their lattice 
sites, V(x;R) will have a more complex landscape than in the T = case, and pc(x|R), 
which is simply an exponential rescaling of ]/(x;R), will again be unable to describe the 



status of the vacancy. The high temperature MD test described in Sec. |III A| provides and 



example of this statement. In Fig. [TT] we compare the classical and quantum densities at the 
configuration of the third panel of Fig.[5| roughly corresponding to the mid- way configuration 
along this local migration event. The classical density was computed at /3 = 0.1, the lowest 
j3 at which its value on crystal minima is negligible. The quantum density is bimodal, with 
one mode localized on each of the crystal sites involved in the migration path. On the hand, 
the classical density is trimodal, with one mode on one crystal site and two modes on the 
other. This test shows that already for this simple system, when the temperature is finite, 
the distribution pc(x|R) is inadequate to represent the state of a vacancy. 
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Appendix B: The case of ab initio force fields. 

Perhaps the most complex extension one can imagine of our representation of a single 
vacancy is to the case of a system interacting via an ab initio force field. In this case, our pro- 
posal is to represent the vacancy in terms of a probe "atom" , with its nucleus and electrons 
equal in charge and number, respectively, to that of the missing atom. Since we assume that 
the Born-Oppenheimer approximation holds, the electrons will not appear explicitly in the 
representation of the vacancy, which will be characterized via the nucleus of the probe atom. 
The electrons of the probe atom will nevertheless contribute to the characterization of the 
vacancy as they determine the effective potential energy of the interaction of the nucleus of 
the probe atom with the nuclei of the real atoms of the sample. Consistently with this, and 
the notation used in the rest of the manuscript, in the following we will denote the positions 
of the nuclei of the atoms and of the probe atom with R and x, respectively. Let us now 
consider the system composed of the nuclei of the physical atoms and of the probe particle, 
and the corresponding electrons in their ground state. Within the Born-Oppenheimer ap- 
proximation the nuclei interact via the potential energy £7(x, R) = (%(r; x, R))# s , where r 
denotes the positions of the real electrons plus the electrons of the probe particle, 'H(r; x, R) 
is the Hamiltonian of the system, and (-) gs denotes the expectation value taken over the 
ground state wavefunction of 7{(r;x, R). We identify the potential V^xjR) of Eq. [TJwith 
£"(x, R). The rest of the vacancy representation is unaffected by the fact that the atoms 
interact via an ab initio force field. 
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